Influence of trapping potentials on the phase diagram of bosonic atoms in optical lattices 



S. M. Giampaolo, F. lUuminati, G. Mazzarella, and S. De Siena 

Dipartimento di Fisica "E. R. Caianiello", Universita di Salerno, INFM UdR di Salerno, 
INFN Sezione di Napoli - Gruppo Collegato di Salerno, 1-84081 Baronissi ( SA ), Italy 

(Dated: November 4, 2004) 

We study the effect of external trapping potentials on the phase diagram of bosonic atoms in optical lattices. 
We introduce a generalized Bose-Hubbard Hamiltonian that includes the structure of the energy levels of the 
trapping potential, and show that these levels are in general populated both at finite and zero temperature. We 
characterize the properties of the superfluid transition for this situation and compare them with those of the 
standard Bose-Hubbard description. We briefly discuss similar behaviors for fermionic systems. 
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I) Introduction. - Exciting progress in the manipulation of 
neutral atoms in optical lattices has recently led, in a series of 
beautiful experiments, to the first verification in atomic sys- 
tems of a quantum ^ase transition from a superfluid to a Mott 
insulator phase [1, 2]. From a more general perspective, the 
physics of ultracold neutral atoms in discrete structures has 
become an ideal testing ground for the study and the realiza- 
tion of complex systems of condensed matter physics l^misl. 

Because of the great current interest in the physics of 
atomic systems in optical lattices, much theoretical effort has 
been devoted to the characterization of the different quantum 
phases arising in these structures. Most of the studies have 
been restricted to situations in which the atoms are assumed 
to be confined in the lowest Bloch band of the periodic lattice 
potential, and to the ground state of the trapping harmonic 
oscillator potential, so that the description is given in terms 
of the standard Bose-Hubbard model However, one 

observes that often, in practical situations, the atoms (either 
fermions or bosons), while remaining confined in the lowest 
band of the periodic lattice potential (for sufficiently low tem- 
peratures), nevertheless may occupy large bundles of excited 
levels of the harmonic confining potential. This is an indica- 
tion that the description in terms of the standard Hubbard and 
Bose-Hubbard models is not always fully adequate, even at 
very low temperatures. 

In the present work we aim to characterize the physics of 
neutral atoms in single-band optical lattices, but including the 
energy structure due to the presence of the trapping poten- 
tial. This leads to the introduction of generalized Hubbard 
and Bose-Hubbard models that are able to explain in a natural 
way the occupation of the excited levels of the trapping po- 
tential. We discuss how such an energy structure affects the 
phase diagram of the system and the nature of the superfluid 
phase. The problem is rather trivial if one considers identiacl 
fermions, but becomes more interesting in the case of bosons. 
We will comment in more detail about fermions in the con- 
clusions. In the following, we consider a dilute ensemble of 
identical, spinless neutral bosonic atoms subject both to a 3-D 
anisotropic harmonic trapping and a 1-D optical lattice po- 
tential. We derive a generalized Bose-Hubbard Hamiltonian 
and we evaluate, both numerically and analytically, the free 
energies of the system. We show that, taking into account 



the possibility to populate the excited levels of the harmonic 
trapping potential, allows the superfluid to become distributed 
among the different energy levels. This effect of distributed 
superfluidity is realized in a single, stable superfluid state, and 
therefore should not be confused with the concept of multi- 
ple (fragmented), unstable superfluid states |7J. This fact al- 
lows to define suitable renormalized quantities that establish 
a precise mapping between the generalized and the standard 
Bose-Hubbard model. 

//) General setting. - The microscopic Hamiltonian for an 
ensemble of bosonic atoms subject to an optical lattice poten- 
tial and confined by an additional, slowly varying, external 
harmonic trapping reads H = T + V + W, where T is the ki- 
netic energy term, V represents the external potential energy 
term and W is the local two-body interaction with coupling 
constant gsB ^ Anh'^aBB/'^m, where m is the mass of the 
atoms and gb b is the s-wave scattering length. In the follow- 
ing, we will always assume boson-boson repulsion a bb > 0. 
The potential energy V is made up by two different contri- 
butions. The first contribution is an harmonic potential Vh 
representing the effects of the quadrupolar trapping magnetic 
field, Vh = mw^ (x^ + A^y^ + A^z^) /2. Here A denotes 
the anisotropy coefficient and uj is the frequency associated to 
the harmonic trap in the x direction. We consider the situation 
A <C 1, representative of the so-called "pancake-shaped" con- 
figuration. The ground state harmonic oscillator length in the 
X direction is = ^hj (jtilo), while for Ly and we have 
that = Ly = Lj./ \f\, i.e. Ly ^, ^ L^. The second contri- 
bution to the potential energy is a 1-D periodic optical lattice 
Vopt — Vb sin^ (ttz/o), where Vb is the maximum amplitude 
of the light shift associated to the intensity of the laser beam 
and a is the lattice spacing related to the wave vector fc of the 
standing laser light by fc = 7r/a. This choice of the geomet- 
rical setting, which is not the one commonly used in exper- 
iments, is merely due to the numerical simplification of the 
energy spectra, and does not affect the general phenomeno- 
logical features, as we will discuss in the following. 

The bosonic field operators can be expanded in the basis of 
the single-particle Wannier wave functions localized at each 
lattice site z,;. The presence of the optical lattice and the 
strong anisotropy define different energy scales along the dif- 
ferent spatial directions. The energy gap between different 
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eigenstates along the y direction is much smaller than the gap 
along the other ones. Because the typical interaction energies 
involved are normally not strong enough in order to excite 
higher vibrational states in the x and z directions, we can re- 
tain only the the lowest vibrational state at each lattice poten- 
tial well. In the harmonic approximation, the Wannier wave 
functions 5'(r) factorize in the product of harmonic oscillator 
states w{r): '^[r) — ^ ai^aw{z — Zi)w°' {y)w{x), where Zi 
is the center of the ith lattice site and 0,^ 0, is the bosonic anni- 
hilation operator acting on the a-th harmonic oscillator level 
at the i-th lattice site. Considering different geometries such 
as a cigar-shaped or a fully isotropic potential will introduce 
extra indices in the orthogonal (x, z) directions, thus compli- 
cating the numerical evaluations but not changing the basic 
physical phenomena. The local ground state spatial extension 

Iz for each lattice potential well is = {a^Eji/ (tt^'Vo)) 
where Eb. — (7r?i)^/2a^m is the lattice recoil energy. 

The condition of anisotropy A ^ 1, implies that one can 
neglect in the Hamiltonian all terms proportional to pow- 
ers of A higher than one. Moreover, in the presence of a 
slowly varying harmonic potential, we may neglect both next- 
to-nearest neighbor hopping and nearest-neighbor interaction 
terms that are some order of magnitude smaller than, respec- 
tively, nearest-neighbor hopping and on-site interaction terms. 
The harmonic Wannier scheme thus leads to the following 
multi-band Bose-Hubbard Hamiltonian 

H = Hi - Y^ta&J,o. , (1) 

<i,j>,a 

Hi = {Eq + Xahw) fli^a +U'Y d'^]5°'ta°'Xl3°'-i;l^i,S • 

Here Hi stands for local terms. We see that neglecting the 
terms proportional to A^ (and higher powers of A) makes the 
Hamiltonian homogeneous. Thus, in this setting (which is 
experimentally realistic), homogeneity is a consequence of 
strong anisotropy. The nonlocal term in Eq. Q is the nearest- 
neighbor hopping contribution with amplitude for a fixed 
harmonic energy level a. In standard experimental situa- 
tions we have Jq > 0, and from the analytical expression 
of the hopping amplitude one finds that Ja+i < J a for every 
a > 0. The first local term in Eq. Q is the sum of local, 
level-dependent energies proportional to the on-site number 
operators „ on the ath harmonic level on site i. This term 
has two contributions: a zero-point energy (i?o) and an excita- 
tion energy {Xahiu) that a boson needs to occupy the ath har- 
monic energy level. The second local term in the Hamiltonian 
is the boson-boson on-site interaction with coupling constant 
U = gBs/ {LxLylz(2Tr)^/'^y The numerical coefficients 

4f readd^f = LyV2^ jr^w'^iy)wP {y)w^ {y)w' iy)dy. 
It is easy to prove that \d"'g \ < 1, with the exception ^q^q — 1. 
Contrary to what happens in the standard Bose-Hubbard set- 
ting 1210001 J we see that the local interaction term involves 
both inter-level interactions and inter-level hoppings which 
do not commute with the on-site number operators n^.Q for 



a given harmonic oscillator level. As a consequence, the two 
local terms of the generalized Bose-Hubbard Hamiltonian Q 
do not in general commute. We then introduce the total on- 
site number operator Ni — fii^a- It is straightforward to 
verify that [Ni, Hi] =0, so that there exists a complete set of 
states that are simultaneous eigenstates of Hi and Ni. There- 
fore, the eigenstates of Hi can be arranged in classes charac- 
terized by a fixed eigenvalue of Ni, with each class possessing 
a lowest eigenvalue. 

Ill) Mean field theory and superfluid structure. - The above 
analysis suggests to divide the problem of finding the eigen- 
states of Hi in different problems in which one determines 
only the eigenstates characterized by a fixed eigenvalue Ni. 
For ease of numerical evaluation it is convenient to adopt 
the grand canonical description K ^ H — fi J^i where 
/i is the chemical potential. We can analyze this problem 
by a mean field approach that corresponds to the approxima- 
tion in which the hopping term is fully decoupled, in anal- 
ogy with the standard Bose-Hubbard problem i^: afad.j,a = 
+ (t>a - 4>i, where 0„ = (a+^) = (a+„) is the 
real (without loss of generality), homogeneous superfluid or- 
der parameter for the ath harmonic energy level. The intro- 
duction of a homogeneous superfluid parameter is justified by 
the homogeneity of the generalized Bose-Hubbard Hamilto- 
nian, which, in the mean-field approximation reads 

Kmf = Hi-n^^Ni 

i 

i,ct i,a 

At fixed values of both the chemical potential and the temper- 
ature, we obtain the eigenvalues of Hi — J^i associated 
to the eigenvectors of Hi. To achieve this goal we proceed 
by truncating at a number rimax of harmonic oscillator lev- 
els, and then diagonalizing the rimax x n-max matrix elements 
of Kmf obtaining the eigenvalues and, hence, the free en- 
ergy. We then iterate the process by increasing Umax- Con- 
vergence is reached when the difference between the iterated 
free energies is lower than a given control parameter This 
is a generalization of the Sheshadri approach to the standard 
Bose-Hubbard problem |8]. A first physical consequence is 
that the mean value of the local occupation number of the a- 
th level is in general nonvanishing, thus explaining the possi- 
ble occupation of the excited levels of the harmonic trapping 
potential, even at zero temperature. In general, the procedure 
outlined above leads to consider several eigenvectors of Hi 
and the large matrices so obtained require thorough numerical 
evaluations. However, with suitable choices of the physical 
parameters only the two lowest eigenstates give non negligible 
contributions, and in this case the analytical diagonalization of 
Kmf is possible. In particular, this situation is realized when 
Jo ^ Xhu! and the chemical potential p is chosen in such a 
way that the state I^q"'*), i.e. the eigenstate with the lowest 
energy in the set of eigenstates of Hi with A'^; = n, and the 
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state I^q"'*'^-'), i.e. the eigenstate with the lowest energy in the 
set of eigenstates of Hi with iVi = n + 1, are nearly degener- 
ate (i.e., the difference between the local eigenvalues of these 
two states is comparable in magnitude with Jq). Let us intro- 
duce the quantity 2 A, the small energy gap between these two 
states. The diagonal terms of the matrix Hamiltonian read 



= -A, 



while the off-diagonal terms due to the hopping are 



(3) 



(4) 



The eigenvalues x 



where Cq. stands for (■0q"-' |act|-0Q"^"'"-' 

of the energy matrix are x 1,2 = iy^A^ + {Y^a Ja(l)aCaf. 
Knowing the eigenvalues, we may write the free energy of 
the system F = In [2 cosh (/3x)] + (Ea -^a-^a) A 

with (3 — (kBT)^^, T the absolute temperature, and fc^ the 
Boltzmann constant. The free energy is a functional of the 
whole set of order (superfluid) parameters (pa determined self- 
consistently through the minimization conditions 

=0 => (j)a^— tanh(/3x) ( Ja4>aCa ) • (5) 

Obviously, this set of self-consistent equations allows for the 
existence of a disordered phase in which all the superfluid or- 
der parameters vanish identically. The crucial question arising 
from Eq.(|5} is whether ordered phases are possible and the 
superfluid can be distributed among the different harmonic 
levels. To have the whole superfluid concentrated on a sin- 
gle level, the set of equations (|5} should allow a solution in 
which only one order parameter, say 4>k with k fixed, is non- 
vanishing. Hence the question is whether a solution of this 
kind is allowed. Let us first consider the case of filling fac- 
tor n = 1, when the states considered are iV'o^'') and iV'o^'')- 
In the state I'i/'o^'')^ due to the presence of exactly one atom 
per lattice site, the boson-boson local interaction vanishes and 
hence |7/)g^^ ) simply factorizes in the product of single-particle 
ground states |l)o, each containing one boson. On the other 
hand, \ipQ ) is not an eigenstate of fii^a for any a, but it may 
be written as a linear combination |1, 1)7,5, characterized by 
having a boson in the 7th level and a boson in the Sth level: 

1^0^^) = J^SK-y ""7,15111 1)7, <5- The coefficients u-y^s are func- 
tions of the relative intensity of the two-body interactions and 

of the gap between the different harmonic levels. Taking into 

(2) 

account the expression of \tpQ ') the coefficients Cq. that ap- 
pear in the off-diagonal terms of the truncated energy matrix 
Eq.© become equal to uo,7 for 7^0 and -\/2 uo,o other- 
wise. From these expressions and from the form of Eqs. (|5j, 
we see that the solution with the superfluid confined in a sin- 
gle band is allowed if and only if the coefficients uq^^ are zero 
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FIG. 1: From top to bottom: the first three nonvanishing superfluid 
order parameters (po, (j)2, and <^4 as a function of the temperature 
kBT/Jo, for Xhuj = 0.1(7 and Ja — Ja+i ~ O.OOIJq. Lines for 
A = 0, points for A — 0.7Jo. Notice the existence of a unique 
critical temperature at fixed values of A. 



for each 7 except at a single fixed value 7 = fc. However, this 



be 



(2) 

condition is inconsistent with the requirement that ' 
an eigenstate of Hi, because this latter constraint implies, for 
every coefficient ^0,7', 



jO,7' 



= (7'Afii^-4'^)"0.7' +4^W7,5C«7 

(5<7 

+ 2\/2^u7,7d°;2;' , (6) 



where 

any 7' such that k 
because if k 
yields 



/,(2) 



For 



is the eigenvalue of Hi associated to \ipo 

7' is even (this assumption is important 
0), solving for uo,fc 



7' is odd, then djj'y 



-1 



Uo.k 



4d 



4 ^7,50?! 



0,k 

0,7' \ 0<<5<7 



rO,7'_ 
7,(5 



2V2^, 

0<7 



a) 



We thus arrive at a set of equations with a number of variables 
less than the number of equations. Because the extra equa- 
tions are linearly independent from each other it is impossible 
to findasetofhamiltonian parameters thatallowswo. 7 = Ofor 
any 7 7^ fc. Numerical evidence shows that qualitatively simi- 
lar results hold for n > 1. In the case n = (hard core limit), 
the superfluid is confined to the harmonic oscillator ground 
state. This result may be explained by observing that in this 
case the local interactions do not play any role. 

IV) Results and discussion. - In Fig. we show the behav- 
ior of the superfluid order parameters associated to different 
harmonic oscillator energy levels as functions of the tempera- 
ture for different values of the energy gap between the states 
IV'o^^) and IV'o^^) (In all fi gures the quantities plotted are di- 
mensionless). From Fig.^we see that the ratio of superfluid 
in the different levels is fixed and independent of the temper- 
ature. This is clear evidence that, although the superfluid is 
distributed among the different levels, the superfluid state is 
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FIG. 2: Critical temperature ksTc/Jo as function of tlie gap A/ Jo 
for different values of Xhu at fixed Jq — Ja+i ~ O.OOlJo. 
Dashed, solid, and dotted lines are respectively for Xhco = 
0.2(7, 0.1(7, 0.01(7. 




FIG. 3: Critical temperature kBTc/Jo as a function of the filling 
factor n for Ja — Ja+i ~ O.OOlJo and Xhui — 0.1(7. The dots cor- 
respond to the analytical solution obtained with the two lowest local 
states. The full line corresponds to the numerical solution including 
the first 40 local states. 



In Fig. |2]we show the behavior of Tc, for varying Xfuv at 
fixed Ja- As far as Tc is concerned, one can see from this 
figure and Eq. (|8} that the transition is analagous to that of 
a standard Bose-Hubbard model with renormalized hopping 
amplitude J_r = ^a^aJa- In Fig.|3lwe show the behavior 
of the critical temperature as a function of the filling factor 
n both for the analytical solution with only the two lowest- 
lying local states and the numerical solution obtained after 
full convergence is reached. The zero-temperature quantum 
phase transition from a superfluid to a Mott-insulator phase is 
recovered for the value Ac at which Tc = 0. From Eq. (|8ji one 
has Ac = Y.a Ja- 



in conclusion, we have studied the properties of an ensem- 
ble of neutral bosonic atoms in an optical lattice, including 
the energy structure due to the presence of a superimposed 
trapping potential. We have shown that the system possesses 
an ordered phase in which the superfluid is distributed among 
the different harmonic energy levels. Our analysis holds in 
a similar way in the much simpler instance of spin-polarized 
fermions and thus provides the physical mechanism explain- 
ing the population of the excited levels of the harmonic trap- 
ping potential. Work is in progress for the case of interact- 
ing fermions, and will be reported elsewhere. An interest- 
ing possibility for future research lies also in the study of 
the energy structure for multicomponent systems and mix- 
tures of bosonic and fermionic atoms in optical lattices. We 
thank INFM, INFN, and MIUR under national project PRIN- 
COFIN 2002 for financial support. 



unique, and therefore stable fl\. From Eq. (|5} the ratios are 
simply: Tq. — Ca/ vSo^' where is the fraction of super- 
fluid in the level a. The fact that the superfluid state is unique 
leads naturally to the existence of a unique critical tempera- 
ture Tc. The latter depends on the energy gap A and on the 
different hopping amplitudes Jq, and may be obtained from 
the set of Eqs. (|5}: it is the temperature for which there exists 
a double solution at (/)ci = for every a: 

/3c A ^ arctanh ( ^ \ , ) . (8) 
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